#include "TMath.h"
#include <iostream>
#include "BinomialErrorInterval.C"
#include "TGraphAsymmErrors.h"
#include "TH1.h"
using namespace std;

//======================================================================
TGraphAsymmErrors * EfficiencyPlot( const TH1 *htot, const TH1 *hsel){
//======================================================================
// -------------------------------------------------------------------------------------------
// Purpose:       Produce TGraph for efficiencyplot with correct Binomial Error interval         
// Author/Date :  Ivo van Vulpen (Ivo.van.Vulpen@nikhef.nl), October 2007                                                                 
// Note:          1) Similar to TGraphAsymmErrors::BayesDivide. Uses private error computation                                       
//                2) Link to BinomialErrorInterval.C                                  
//                   http://www.nikhef.nl/~ivov/BinomialErrorInterval/BinomialErrorInterval.C    
// -------------------------------------------------------------------------------------------
 
  Int_t Nbins_tmp = (Int_t)htot->GetNbinsX(); 
  const Int_t Nbins = Nbins_tmp;

  // Fill Arrays
  Double_t X[Nbins];    Double_t X_el[Nbins];    Double_t X_eh[Nbins];
  Double_t Eff[Nbins];  Double_t Eff_el[Nbins];  Double_t Eff_eh[Nbins];

  Double_t Ntot;
  Double_t Nsel;
  Double_t ProbRegion[2];

  Int_t id = 0;

  // Loop over all bins and compute efficiency and it's uncertainty
  for(Int_t i_bin = 1; i_bin< Nbins+1; i_bin++){

    // Compute Efficiencies
    Ntot  = (Double_t)htot->GetBinContent(i_bin);
    Nsel  = (Double_t)hsel->GetBinContent(i_bin);

    // skip bins w/o content
    if ( Ntot > 0 ) {

        // Compute X-values
        X[id]    = htot->GetBinCenter(i_bin); 
        X_el[id] = 0.;
        X_eh[id] = 0.;

        cout << Nsel << ":" << Ntot << endl;

        BinomialErrorInterval(Nsel,Ntot,ProbRegion);
        Eff[id]    = ( Ntot != 0) ?  (Double_t)Nsel/(Double_t)Ntot : 0.;
        Eff_el[id] = ( Eff[id] - ProbRegion[0] );
        Eff_eh[id] = ( ProbRegion[1] - Eff[id] );
 
        // Correct strange things
        if( Nsel < 0. || Nsel > Ntot || Ntot < 0.){
          cout << " Strange bin: X = " << X[id] << "   Nsel = " << Nsel << "   Ntot = " << Ntot << "  --> setting efficiency to 0." << endl;
          Eff[id]    = 0.;
          Eff_el[id] = 0.;
          Eff_eh[id] = 0.;
        }

        id++;
    }
  } // end loop over bins


  // Set some specifics for the Graph and return it to user
  TGraphAsymmErrors *graph = new TGraphAsymmErrors(id,X,Eff,X_el,X_eh,Eff_el,Eff_eh);
  graph->SetTitle("Efficiency Plot");
  graph->SetMarkerColor(4);
  graph->SetMarkerStyle(20);

  return graph;

} // end EfficiencyPlot

